NAVAL 

POSTGRADUATE 

SCHOOL 

MONTEREY,  CALIFORNIA 


THESIS 


ADAPTIVE  SELECTIONS  OF  SAMPLE  SIZE  AND 
SOLVER  ITERATIONS  IN  STOCHASTIC  OPTIMIZATION 
WITH  APPLICATION  TO  NONLINEAR  COMMODITY 
FLOW  PROBLEMS 

by 

David  A.  Vondrak 

March  2009 

Thesis  Advisor: 
Second  Reader: 

Johannes  0.  Roy  set 

R.  Kevin  Wood 

Approved  for  public  release;  distribution  is  unlimited 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  No.  0704-0188 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instruction, 
searching  existing  data  sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send 
comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing  this  burden,  to 
Washington  headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA 
22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188)  Washington  DC  20503. _ 

2.  REPORT  DATE 

March  2009 


6.  AUTHOR(S)  David  A.  Vondrak 


11.  SUPPLEMENTARY  NOTES  The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official  policy 
or  position  of  the  Department  of  Defense  or  the  U.S.  Government. 


13.  ABSTRACT  (maximum  200  words) 

We  present  an  algorithm  to  approximately  solve  certain  stochastic  nonlinear  programs  through  sample- 
average  approximations.  The  sample  sizes  in  these  approximations  are  selected  by  approximately  solving  optimal- 
control  problems  defined  on  a  discrete-time  dynamic  system.  The  optimal-control  problem  seeks  to  minimize  the 
computational  effort  required  to  reach  a  near-optimal  objective  value  of  the  stochastic  nonlinear  program.  Unknown 
control-problem  parameters  such  as  rate  of  convergence,  computational  effort  per  solver  iteration,  and  optimal  value 
of  the  program  are  estimated  within  a  receding  horizon  framework  as  the  algorithm  progresses.  The  algorithm  is 
illustrated  with  single-commodity  and  multi-commodity  network  flow  models.  Measured  against  the  best  alternative 
heuristic  policy  we  consider  for  selecting  sample  sizes,  the  algorithm  finds  a  near-optimal  objective  value  on  average 
up  to  17%  faster.  Further,  the  optimal-control  problem  also  leads  to  a  40%  reduction  in  standard  deviation  of 
computing  times  over  a  set  of  independent  runs  of  the  algorithm  on  identical  problem  instances.  When  we  modify  the 
algorithm  by  selecting  a  policy  heuristically  in  the  first  stage  (only),  we  improve  computing  time,  on  average,  by 
nearly  47%  against  the  best  heuristic  policy  considered,  while  reducing  the  standard  deviation  across  the  independent 
runs  by  55%. 


16.  PRICE  CODE 


NSN  7540-0 1  -280-5500  Standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Std.  239-18 


20.  LIMITATION  OF 
ABSTRACT 


15.  NUMBER  OF 
PAGES 

57 


14.  SUBJECT  TERMS  Nonlinear  Stochastic  Optimization,  Optimal  Control,  Dynamic 
Programming,  Network  Commodity  Flow,  Sample  Average  Approximation,  Projected 
Gradient  Method 

18.  SECURITY 
CLASSIFICATION  OF  THIS 
PAGE 

Unclassified 


19.  SECURITY 
CLASSIFICATION  OF 
ABSTRACT 

Unclassified 


17.  SECURITY 
CLASSIFICATION  OF 
REPORT 

Unclassified 


12b.  DISTRIBUTION  CODE 


12a.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 
Approved  for  public  release;  distribution  is  unlimited 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Naval  Postgraduate  School 
Monterey,  CA  93943-5000 

9.  SPONSORING  /MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

N/A 


5.  FUNDING  NUMBERS 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 


4.  TITLE  AND  SUBTITLE  Adaptive  Selections  of  Sample  Size  and  Solver 
Iterations  in  Stochastic  Optimization  with  Application  to  Nonlinear  Commodity 
Flow  Problems 


3.  REPORT  TYPE  AND  DATES  COVERED 

Master’s  Thesis 


1.  AGENCY  USE  ONLY  (Leave  blank) 


1 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


11 


Approved  for  public  release;  distribution  is  unlimited 


ADAPTIVE  SELECTIONS  OF  SAMPLE  SIZE  AND  SOLVER  ITERATIONS  IN 
STOCHASTIC  OPTIMIZATION  WITH  APPLICATION  TO  NONLINEAR 
COMMODITY  FLOW  PROBLEMS 


David  A.  Vondrak 

Lieutenant  Commander,  United  States  Navy 
B.S.,  United  States  Naval  Academy,  1997 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


MASTER  OF  SCIENCE  IN  OPERATIONS  RESEARCH 


from  the 


NAVAL  POSTGRADUATE  SCHOOL 
March  2009 


Author:  David  A.  Vondrak 


Approved  by:  Johannes  O.  Royset 

Thesis  Advisor 


R.  Kevin  Wood 
Second  Reader 


Robert  F.  Dell 

Chairman,  Department  of  Operations  Research 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


IV 


ABSTRACT 


We  present  an  algorithm  to  approximately  solve  certain  stochastic  nonlinear 
programs  through  sample-average  approximations.  The  sample  sizes  in  these 
approximations  are  selected  by  approximately  solving  optimal-control  problems  defined 
on  a  discrete -time  dynamic  system.  The  optimal-control  problem  seeks  to  minimize  the 
computational  effort  required  to  reach  a  near-optimal  objective  value  of  the  stochastic 
nonlinear  program.  Unknown  control-problem  parameters  such  as  rate  of  convergence, 
computational  effort  per  solver  iteration,  and  optimal  value  of  the  program  are  estimated 
within  a  receding  horizon  framework  as  the  algorithm  progresses.  The  algorithm  is 
illustrated  with  single-commodity  and  multi-commodity  network  flow  models.  Measured 
against  the  best  alternative  heuristic  policy  we  consider  for  selecting  sample  sizes,  the 
algorithm  finds  a  near-optimal  objective  value  on  average  up  to  17%  faster.  Further,  the 
optimal-control  problem  also  leads  to  a  40%  reduction  in  standard  deviation  of 
computing  times  over  a  set  of  independent  runs  of  the  algorithm  on  identical  problem 
instances.  When  we  modify  the  algorithm  by  selecting  a  policy  heuristically  in  the  first 
stage  (only),  we  improve  computing  time,  on  average,  by  nearly  47%  against  the  best 
heuristic  policy  considered,  while  reducing  the  standard  deviation  across  the  independent 
runs  by  55%. 


v 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


vi 


TABLE  OF  CONTENTS 


I.  INTRODUCTION . 1 

A.  BACKGROUND . 1 

B.  LITERATURE  REVIEW . 2 

C.  RESEARCH  GOAL . 3 

D.  STRUCTURE  OF  THESIS  AND  CHAPTER  OUTLINE . 4 

II.  PROBLEM  DEFINITION  AND  FORMULATION . 5 

A.  PROBLEM  FORMULATION . 5 

1.  Expected-value  Problem . 5 

2.  Approximating  Problem . 6 

3.  Properties  of  Sample  Average . 7 

B.  ALGORITHM . 8 

C.  PRECISION-ADJUSTMENT  PROBLEM . 9 

III.  METHODOLOGY . 11 

A.  PRECISION  ADJUSTMENT  CONTROL . 11 

1.  Controlling  Sample  Size . 11 

2.  Dynamic  Program . 14 

B.  PARAMETER  ESTIMATION . 17 

1 .  Estimating  V ariance . 18 

2.  Estimation  of  Rate-of-Convergence  Coefficient  and  Optimal 

Objective  Value . 18 

C.  STOPPING  CRITERION . 20 

IV.  COMPUTATIONAL  STUDY . 23 

A.  SINGLE-COMMODITY  NETWORK  FLOW . 23 

B.  MULTI-COMMODITY  NETWORK  FLOW . 25 

C.  COMPUTATIONAL  STUDY . 27 

V.  CONCLUSIONS . 35 

LIST  OF  REFERENCES . 37 

INITIAL  DISTRIBUTION  LIST . 39 


vii 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


LIST  OF  FIGURES 


Figure  1 .  Expected  Value  function  compared  to  Sample  Average  function . 6 

Figure  2.  Approximate  Normal  Distribution  of  fN(x) . 8 

Figure  3.  Approximate  Normal  Distribution  of  fN  (jc*  ) . 12 

Figure  4.  Approximate  Truncated  Normal  Distribution  for  / (x^  )  . 13 

Figure  5.  Truncated  normal  distributions  of  function  values . 14 

Figure  6.  Representative  illustration  of  dynamic  program . 1 5 

Figure  7.  Discretization  of  truncated  nonnal  with  small  Nk  and  nk . 17 

Figure  8.  Discretization  of  truncated  nonnal  with  large  Nk  and  nk  . 17 

Figure  9.  Transportation  Grid  Network . 24 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


x 


LIST  OF  TABLES 


Table  1.  Average  and  standard  deviation  of  computing  times  of  CA  with  Model- 

Predictive  Control  (MPC)  policies  and  heuristic  policies  applied  to  SCF . 30 

Table  2.  Average  and  standard  deviation  of  computing  times  of  CA  with  Model- 

Predictive  Control  (MPC)  policies  and  heuristic  policies  applied  to  MCF . 32 


xi 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


EXECUTIVE  SUMMARY 


Optimization  of  stochastic  programs  is  challenging  in  part  because  there  is  no 
closed-form  solution,  and  because  solution  algorithms  tend  to  be  computationally 
intensive.  The  difficulty  arises  because  the  objective  function  and/or  constraint 
functions  are  given  by  expectations  that  cannot  be  evaluated  exactly.  Hence, 
approximations  are  usually  employed  in  optimization  algorithms  to  estimate  such 
functions.  One  class  of  such  algorithms  uses  sample-average  approximations  within  a 
nonlinear  approximating  problem  (AP).  However,  there  is  no  straightforward  means  to 
determine  a  sample  size  for  use  in  these  algorithms.  Most  algorithms  resort  to  using 
heuristic  policies  for  determining  an  appropriate  sample  size.  In  addition  to  sample  size, 
because  the  AP  is  nonlinear,  a  suitable  number  of  iterations  of  a  nonlinear  programming 
solver  must  also  be  selected  for  solving  it.  Most  often,  the  number  of  solver  iterations  is 
selected  prior  to  calculations  beginning.  It  is  often  difficult  in  practice  to  select  sample 
sizes  and  number  of  solver  iterations  that  balances  accuracy  and  computational  effort. 

We  develop  a  discrete-time  dynamic  system,  the  optimal  control  of  which 
determines  a  policy  for  sample  size  and  a  number  of  solver  iterations  for  use  in  the  AP. 
The  optimal-control  problem  seeks  to  minimize  the  computational  effort  required  to 
reach  a  near-optimal  objective  value  of  the  stochastic  nonlinear  program.  The  optimal- 
control  problem  is  approximately  solved  within  a  receding  horizon  framework,  allowing 
repeated  estimation  of  unknown  parameters. 

Empirical  studies  are  performed  on  nonlinear  single  and  multiple-commodity 
network-flow  problems  on  a  grid  constructed  of  50  nodes  with  134  arcs.  The  optimal- 
control  problem  selects  sample  sizes  and  solver  iterations  that  lead  to  near-optimal 
objective  values  in  less  time  than  alternative  heuristic  policies  in  all  instances  tested. 
Measured  against  the  best  alternative  policy  we  consider  for  selecting  sample  sizes,  the 
algorithm  finds  a  near-optimal  objective  value  on  average  up  to  17%  faster.  Further,  the 
optimal-control  problem  also  leads  to  a  40%  reduction  in  standard  deviation  of 
computing  times  over  a  set  of  independent  runs  of  the  algorithm  on  identical  problem 


instances.  The  unknown  parameters  in  the  optimal-control  problem  may  be  poorly 
estimated  prior  to  the  first  stage  of  the  algorithm,  which  may  result  in  a  poor  policy  for 
the  first  stage.  When  we  modify  the  algorithm  by  selecting  a  policy  heuristically  in  the 
first  stage  (only),  we  improve  computing  time,  on  average,  by  nearly  47%  against  the 
best  heuristic  policy  considered,  while  reducing  the  standard  deviation  across  the 
independent  runs  by  more  than  one-half. 
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I.  INTRODUCTION 


A.  BACKGROUND 

Optimization  of  stochastic  programs  has  become  a  focus  of  much  research  over 
the  past  decade.  These  problems  are  of  particular  interest  in  part  because  they  have  no 
closed-form  solutions  and  solution  algorithms  tend  to  be  computationally  intensive.  The 
difficulty  arises  because  the  objective  function  and/or  constraint  functions  are  given  by 
expectations  that  cannot  be  evaluated  exactly.  Hence,  approximations  are  usually 
employed  to  estimate  such  functions.  Approximations  may  provide  precise,  high-fidelity 
estimates  with  high  computational  cost,  or  may  provide  low  fidelity  with  less 
computational  cost.  This  research  focuses  on  finding  a  balance  between  these  two 
important  factors. 

Numerous  design  and  planning  applications  require  the  optimization  of 
stochastic-programming  problems.  In  aerodynamics,  various  problems  arise  in  the 
optimization  of  a  three-dimensional  wing  design  (Alexandrov  et  ah,  2001).  In  one  civil- 
engineering  discipline,  problems  arise  from  structural  optimization  of  bridges  or  support 
structures  subject  to  failure  probability  constraints  (Polak  and  Royset,  2007;  Royset  and 
Polak,  2007).  Such  problems  may  seek  to  optimize  the  cross-sectional  dimensions  of  a 
support  column  consisting  of  a  material  of  particular  yield  strength  subjected  to  various 
bending  moments.  Structural  loading  of  the  support  column  weighs  heavily  on  design 
considerations  due  to  inherent  failure  probabilities  of  the  materials.  Many  additional 
examples  of  stochastic-programming  problems  can  be  found  in  the  areas  of  logistics  and 
supply-chain  planning.  Numerous  papers  discuss  solution  methodologies  for  optimal 
design  of  supply-chain  management  including,  but  not  limited  to  productions  lines, 
consumer  demand,  availability  of  raw  materials,  warehousing  and  transportation  of 
goods.  In  Santoso  et  al.  (2003),  a  proposed  stochastic-programming  model  and  solution 
algorithm  are  used  to  compute  high-quality  solutions  to  large-scale  stochastic  supply 
chain  design  problems.  Poojari  et  al.  (2006)  presents  a  method  to  solve  discrete  resource 
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allocation  problems  in  the  presence  of  future  uncertainties  for  supply-chain  planning  and 
Sox  and  Muckstadt  (1997)  describe  a  finite-horizon  stochastic  optimization  model  for  a 
stochastic  lot-scheduling  problem. 

One  approach  to  approximately  solve  optimization  problems  defined  in  terms  of 
expectations  is  to  use  Sample  Average  Approximations  (SAA),  e.g.,  (Ruszczynski  and 
Shapiro,  2007).  That  is,  one  or  more  approximating  problems  (APs)  are  solved, 
problems  that  replace  each  expectation  with  a  standard  sample  average  approximation. 
The  AP  may  be  viewed  as  a  deterministic  mathematical  program,  a  nonlinear 
mathematical  program  in  our  case.  A  very  large  sample  size  in  the  sample-average 
approximation  will  provide  a  precise  estimate  of  the  expected-value  function,  but  the 
computational  effort  required  to  solve  the  AP  for  this  sample  size  may  be  far  too  high. 
As  an  alternative,  this  research  examines  the  computational  effort  associated  with  solving 
a  sequence  of  APs  where  an  efficient  sample  size  is  chosen  at  each  stage  of  the  sequence. 
Coarse  approximations  are  made  early  and,  as  the  calculations  evolve,  adaptive 
adjustments  to  the  sample  size  are  made,  increasing  the  precision  of  the  results.  This 
adaptive  control  strategy  is  intuitively  appealing  as  gains  towards  optimality  come 
initially  at  low  computational  cost  through  coarse  approximations,  and  fidelity  is 
increased  as  larger  samples  are  used  near  the  completion  of  the  algorithm. 

In  addition  to  finding  an  efficient  sample  size  for  use  with  calculations,  we  are 
also  concerned  with  selecting  an  efficient  number  of  iterations  for  the  chosen  nonlinear 
programming  (NLP)  solver  used  in  solving  various  instance  of  the  AP.  The  number  of 
solver  iterations  has  a  profound  impact  on  computational  cost  as  well.  Our  approach  will 
determine  an  efficient  number  of  iterations  for  the  NLP  solver  to  carry  out  along  with  an 
efficient  sample  sizes  for  defining  the  APs  as  the  overall  algorithm  progresses. 

B.  LITERATURE  REVIEW 

Stochastic  programming  solution  methodologies  incorporate  both  mathematical- 
programming  techniques  and  statistically  motivated  approximations.  These 
approximations  are  generally  internal  or  external  in  nature.  The  distinction  is  found  in 
the  management  of  sample  selection.  With  external  methodologies  (Royset  and  Polak, 

2 


2007),  samples  are  taken  before  the  solution  algorithm  begins  (or  they  could  be  taken 
before  the  algorithm  begins),  and  no  additional  sampling  is  performed  during  the 
optimization  process.  In  internal  methods  (Higle  and  Zhao,  2004),  samples  are  intrinsic 
to  the  iterative  solution  process,  perfonned  whenever  the  algorithm  requires  the 
estimation  of  expectations.  An  alternative  to  an  external  sampling  approach  with  a  fixed 
sample  size  is  to  vary  the  sample  size  during  the  calculations  using  closed-loop  or  open- 
loop  techniques  (Polak  and  Royset,  2007).  With  a  closed-loop  technique,  the  sample  size 
is  adapted  during  the  calculations.  For  example,  if  the  objective  value  in  a  given  iteration 
falls  below  a  moving  floor,  the  sample  size  may  be  increased.  On  the  other  hand,  in  an 
open-loop  technique,  sample  sizes  are  preassigned  (Polak  and  Royset,  2007).  In  most 
cases,  the  sample  sizes  increase  as  iterations  progress  towards  an  optimal  solution. 

He  and  Polak  (1990)  describe  a  method  for  handling  progressively  finer  stages  of 
discretization  for  semi-infinite  optimization  problems,  which  uses  a  precision-adjustment 
strategy  relevant  to  the  present  work.  They  fonnulate  an  auxiliary  optimization  problem 
that  determines  the  number  of  solver  iterations  and  precision  level  at  different  stages  of 
the  calculations  so  that  overall  computing  time  is  minimized  approximately.  The  present 
work  is  motivated  by  this  study. 

C.  RESEARCH  GOAL 

The  goal  of  this  research  is  to  find  an  efficient  way  of  selecting  sample  size  and 
number  of  solver  iterations  for  use  in  solving  a  nonlinear  stochastic  program  through  the 
solution  of  a  sequence  of  APs  that  use  sample  average  approximations.  We  formulate  a 
discrete-time  optimal-control  problem  that  we  solve  approximately  to  obtain  a  precision- 
adjustment  policy  for  detennining  the  sample  size  and  number  of  solver  iterations.  This 
policy  makes  coarse  approximations  in  the  early  stages  of  the  problem  and,  as 
progression  to  the  optimal  objective  value  is  achieved,  the  sample  size  increases  to  ensure 
convergence  to  a  locally  optimal  solution.  The  optimal-control  problem  is  formulated 
with  the  objective  to  minimize  the  amount  of  computational  work  necessary  to  reach  a 
locally  optimal  solution  to  the  AP. 
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D. 


STRUCTURE  OF  THESIS  AND  CHAPTER  OUTLINE 


This  thesis  is  organized  into  five  chapters  including  the  Introduction.  Chapter  II 
provides  a  discussion  on  the  formulation  of  a  nonlinear  stochastic  problem,  an  outline  of 
the  algorithm  we  intend  to  use  to  solve  stochastic  nonlinear  network  flow  problems  and 
introduces  our  precision-adjustment  problem.  Chapter  III  introduces  the  methodology  we 
propose  to  use  for  selecting  an  efficient  sample  size  and  number  of  solver  iterations  in  the 
solution  algorithm;  stopping  criterion  are  also  discussed.  Chapter  IV  formulates  and 
solves  two  types  of  stochastic  nonlinear  network  flow  problems.  It  compares  our 
algorithm’s  efficiency  with  the  efficiency  of  similar  algorithms  that  use  heuristic  sample- 
size  and  solver-iteration  policies.  Chapter  V  summarizes  the  research  and  presents  the 
main  findings  and  insights. 
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II.  PROBLEM  DEFINITION  AND  FORMULATION 


This  chapter  formulates  a  stochastic  nonlinear  program  and  presents  a  conceptual 
algorithm  we  use  as  our  solution  approach.  Further,  it  describes  the  precision-adjustment 
problem  we  face  within  the  conceptual  algorithm. 

A.  PROBLEM  FORMULATION 

This  section  describes  the  formulation  of  the  “expected-value  problem”  and  an 
approximating  problem  which  will  be  used  as  a  surrogate  to  approximately  solve  the 
expected-value  problem. 

1.  Expected-value  Problem 

Consider  the  optimization  of,  say,  an  engineering  design  or  logistics  problem 
defined  by  the  random  function  F(x,  co) ,  where  x  e  X  cz  \  "  is  a  vector  of  continuous 
decision  variables  and  co  is  a  vector  of  continuous  random  variables  defined  on  the 
probability  space  (Q ,F,P) .  This  situation  results  in  a  difficult  stochastic  optimization 
problem  of  the  form: 

EP  :  min  {/ (jc)  :=  E[F(x,  co)] } ,  (2.1) 

xeAczj  n 

where  X  is  a  compact  subset  of  ;  ”  representing  feasible  solutions  of  the  problem,  and  E 
is  the  expectation  taken  with  respect  to  the  known  probability  distribution  P  of  the 
random  vector  co .  We  assume  that,  for  every x  e  X  ,  the  expected-value  function  is  well 
defined  and  smooth  (continuously  differentiable).  Further,  we  denote  the  optimal  value 
of  (2.1)  as  /*  with  a  set  of  optimal  solutions  denoted  byX*.  Further  we  indicate  a  set  of 
^-optimal  solutions  by  X*,  i.e.,  for  any  e  >  0 

X;:={xeX\f(x)-r<e}.  (2.2) 

Because  the  distribution  of  F{x,co)  is  unknown,  we  cannot  compute  the 
expectation  in  closed  form,  so  we  approximate  the  expected-value  function  by  the 
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sample-average  estimator  discussed  below.  Uncertainty  may  also  be  introduced  in  the 
constraint  functions  defining  X;  however,  our  research  focuses  only  on  uncertainty  in  the 
objective  function. 

2.  Approximating  Problem 

To  approximate  f{x)  =  EfFXv,  a>)]  we  use  a  sample  average  defined  by: 

(2-3) 

A  j= i 

where  co1  ,co2,...,coN ,  is  a  sample  of  size  N  consisting  of  independent,  identically 
distributed  (iid)  random  variables.  Moreover,  we  define  an  approximating  problem  (AP): 

AP:  min  fN(x) ,  (2.4) 

xeX 

with  an  optimal  value  denoted  by  f*.  We  refer  to  (2.1)  and  (2.4)  as  the  “expected- 
value”  and  “approximating  problems,”  respectively. 

The  expected-value  function,  for  instance,  may  be  represented  by  the  function 
depicted  by  the  solid  line  in  Figure  1.  This  function,  as  stated  above,  cannot  be  computed 
in  closed  form,  so  we  approximate  the  function’s  form  by  the  sample  average  as  shown  in 
Figure  1.  As  discussed  in  Section  3  below,  we  expect  the  sample  average  to  take  the 
form  of  the  expected-value  function  as  the  sample  size  approaches  infinity,  ft  is  known 
that  solving  the  AP  with  an  appropriate  sample  size  provides  a  reasonable  approximation 
to  the  solution  to  EP  (Ruszczynski  and  Shapiro,  2007). 


Figure  1.  Expected  Value  function  compared  to  Sample  Average  function 
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3.  Properties  of  Sample  Average 


Monte  Carlo  simulation  can  be  used  to  generate  an  iid  sample  of  oj's  of  size  N  to 
obtain  F(x,col),F(x,co1),...,F(x,coN)  for  use  in  (2.3).  It  is  well  known  that  fN(x)  is  an 
unbiased  estimator  of  f(x)  i.e., 


E  [fN(x)]  =  f(x), 

Var [fN(x)\  =  a\x)/  N, 

where 

cr2  (x)  =  Var(F(x,  oo))  . 

Because  the  generated  sample  is  iid  and  given  rather  weak  general  assumptions,  it 
follows  from  the  Law  of  Large  Numbers  (LLN)  that  fN(x)  converges  pointwise  to  fix) 
with  probability  one,  as  N  — >  oo  (Ruszczynski  and  Shapiro,  2007).  Therefore,  we  would 
expect  the  optimal  value  and  optimal  solution  of  the  AP  to  converge  to  those  of  EP,  as 
N  — >  oo  .  This  is  made  precise  in  the  following  proposition. 

Proposition  1.  (Prop.  4.2;  Ruszczynski  and  Shapiro,  2007).  If  the  pointwise 
LLN  holds,  i.e.,  fN(x)  converges  to  fix)  uniformly  on  X,  with  probability  one 

as  N  — »  oo  ,  then  f*  converges  to /*  with  probability  one,  as  N  — >  oo  . 

Moreover,  if  E[F(x, co)2~\<  oo ,  then  it  follows  under  weak  assumption  from  the 
central  limit  theorem  (CLT)  that  4n (fN(x)~  f  (x))  =>  Yx,  where  =>  denotes 
convergence  in  distribution  and  Yx  is  a  zero-mean  normal  random  variable  with  variance 
<j2(x)  (Ruszczynski  and  Shapiro,  2007).  Hence,  for  a  large  A,  fN{x)  is  approximately 
normally  distributed  with  mean  /  (x),  and  variance  cr2(x)  /  N  .  Figure  2  illustrates  this 
result. 

We  know  that 


min 

xeX 


E[/„(*)]>E 


"  "  "  /  ',  !■  >■> 

xeX 
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and,  because  fN{x)  is  an  unbiased  estimator  of  /(x)  ,  it  follows  that  E [f* ]  < /* .  That 

is,  the  AP’s  optimal  objective  value  is  a  downward-biased  estimate  of  the  optimal 
objective  value  of  the  EP. 


Figure  2.  Approximate  Normal  Distribution  of  fN(x) 


B.  ALGORITHM 

We  consider  a  conceptual  algorithm  of  the  following  form  for  approximately 
solving  EP. 

Conceptual  Algorithm  (CA) 

Data:  Optimality  tolerance  e  >  0;  initial  point  ijel, 

Step  1 :  Set  stage  counter  k=  1 . 

Step  2:  Determine  Nk  and  nk. 

Step  3:  Carry  out  nk  iterations  of  a  solver  applied  to: 

min  fN  (x) ,  started  with  xk0  using  the  N,  found  at  Step  2. 

xeX  k 

Step  4:  If  x*  e  X*  ,  then  Stop.  Else  set  x{'t '  =  x*  ,  replace  k  by  k+l  and 
go  to  Step  2. 
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We  refer  to  Steps  2-4  as  a  “stage.”  CA  uses  a  single  set  of  sample  realizations  for 
each  stage,  but  all  samples  are  independent  between  stages.  For  purposes  of  this 
research,  we  have  chosen  the  Projected  Gradient  Method  (PGM)  as  our  nonlinear 
programming  solver  (for  example,  see  Polak  1997,  p.  66).  Flowever,  it  is  worth 
mentioning  that  any  linearly  convergent  nonlinear  programming  solver  may  be  used.  Our 
goal  is  to  determine  Nk  and  nk  for  each  stage  that  approximately  minimizes  the  total 
computational  effort  required  by  CA  to  reach  a  near-optimal  objective  value  to  EP.  We 
refer  to  the  rule  for  selecting  Nk  and  nk  as  a  “policy.” 

C.  PRECISION-ADJUSTMENT  PROBLEM 

We  define  the  task  of  selecting  the  appropriate  Nk  and  nk  for  each  stage  of  the 

CA  as  the  Precision- Adjustment  Problem  (PAP).  We  formulate  PAP  as  a  particular 
discrete-time  optimal-control  dynamic  program,  which  is  discussed  in  the  next  chapter. 
PAP  is  solved  approximately  using  dynamic  programming  to  obtain  a  precision- 
adjustment  policy  for  detennining  Nk  and  n, .  This  policy  adaptively  adjusts  the  sample 

size  and  number  of  solver  iterations  between  stages.  PAP  is  fonnulated  to  minimize  the 
amount  of  computational  work  necessary  to  reach  a  near-optimal  objective  value  of  the 
stochastic  nonlinear  program.  After  carrying  out  nk  solver  iterations,  we  abandon  further 

progress  to  the  locally  optimal  solution  for  the  current  AP,  and  use  this  iterate  as  a  wann 
start  for  the  next  stage  of  CA. 
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III.  METHODOLOGY 


This  chapter  describes  how  we  control  the  precision  within  the  conceptual 
algorithm  and  how  we  estimate  the  parameters  for  the  precision  adjustment.  Further,  the 
chapter  presents  the  criterion  by  which  we  stop  the  conceptual  algorithm  (in  Step  4). 

A.  PRECISION  ADJUSTMENT  CONTROL 

This  section  presents  a  discrete-time  dynamic  system,  the  optimal  control  of 
which  will  determine  a  sample  size  Nk  and  a  number  of  solver  iterations  nk  to  be  used 
for  each  stage  of  CA. 

1.  Controlling  Sample  Size 

The  dynamic  system  described  in  Royset  (2009)  lays  the  foundation  for  this 
research.  Following  his  approach,  we  first  identify  the  asymptotic  distributions  of  the 
progress  made  by  CA.  We  assume  that  the  solver  used  in  Step  3  of  CA  is  unifonnly 
linearly  convergent  (see  Royset,  2009),  i.e.,  there  exists  a  6*e(0,l)  such  that 

fN {PN (*)) - /;  <  0{fN (x) - /; )  for  all  rel  and  Ae¥,  where  ¥=  (1,2,3, ...}and 
PN(x)  is  the  iterate  found  after  carrying  out  one  iteration  of  the  solver  used  in  the  CA 
starting  from  x  with  sample  size  N .  We  refer  to  9  as  the  rate-of-convergence 
coefficient.  We  let  P"1  (x)  denote  the  iterate  from  the  solver  found  at  stage  k,  after  nk 

iterations  with  a  sample  of  size  Nk  starting  from  x.  It  follows  from  the  assumption  of 
linear  convergence  and  from  the  optimality  of  f*N  ,  that  for  any  xel, 

/;  *  fNk  (*»  *  rNk + °nk  uNt  w  -  /; ) 

with  probability  one.  Moreover,  based  on  rather  weak  assumptions  and  theorems 
introduced  in  Royset  (2009),  if  Nk  and nk  are  large,  fNk(xk„k)  is  approximately 
distributed  as 
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•W*  +  0'h  (/(<;! )-n,cy\x)INk),  (3.1) 

where  Af(/j,cr)  represents  a  normally  distributed  random  variable  with  mean  //  and 
variance  cr .  Figure  3  illustrates  a  possible  approximate  distribution  of  fN  (x*  )  as  stated 
from  this  conclusion. 


Further,  Royset  (2009)  shows  that  we  can  heuristically  approximate  the 
distribution  /(jc*  )  with  truncation  at  /*  to  account  for  the  fact  that 

/ (x)  >  /*  for  all  x  e  X  .  Therefore  our  approximate  distribution  for  f(P„  (jc*_1  )) 
conditional  on  xk~x  is: 

nk- 1 

w* + (/(<:!  )-rW{x)iNk,n,  o.d 

where  2V  denotes  a  truncated  normally  distributed  random  variable  based  on 
JV(f*  +  0"k (f  (xkn~\ )-/*), cr2(x*)/ Nk),  with  a  truncation  threshold  /* .  Figure  4 
illustrates  this  distribution.  As  discussed  in  Subsection  2,  we  can  control  this  distribution 
by  changing  Nk  and  nk . 
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Figure  4.  Approximate  Truncated  Normal  Distribution  for  / (jc*  ) 

Since  (3.2)  approximately  holds  for  any  k,  we  can  use  that  expression  to  estimate 
the  quality  of  a  solution  generated  by  the  CA  after  any  number  of  stages  given  a  selection 
of  sample  sizes  and  number  of  iterations.  The  situation  is  illustrated  in  Figure  5,  where 
the  f(x°n)  and  a  selection  of  (TV,,  «,)  determines  the  probability  distribution  of  f(x)h ) . 

Given  an  outcome  of  that  distribution,  a  selection  of  (N2,  n2 ) ,  gives  the  distribution  of 
/(*!,),  etc. 

The  values  for  /*,  9,  and  a2(x  )  in  the  distribution  (3.2)  are  unknown,  however. 
Since  x*t  |  is  known  at  the  beginning  of  the  kth  stage,  we  can  estimate  /(x*”j),  and  we 
construct  estimation  schemes  as  shown  in  Subsection  B  below  to 
estimate  /*,  9 ,  and  a2(x' ) .  We  now  develop  a  dynamic  program  as  discussed  in  Royset 
(2009). 
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f(x\h)  is  a  truncated  normal  with  distribution 

a* +0-(/(<)- n,°\x)i  nx,d 


f  (xf  )  is  a  truncated  normal  with  distribution 

(/* n2,d 

f  (xl  )  is  a  truncated  normal  with  distribution 

_ (r+d'Hnxi)-rW(x)/N3,n _ 

Figure  5.  Truncated  normal  distributions  of  function  values 

2.  Dynamic  Program 

Beginning  at  stage  k,  we  define  estimates  of  f(xk~l),f*,  6 ,  and  cr2(x*)  as 
rf,  r* ,  rg,  and  r  2  respectively  and  as  in  Royset  (2009),  use  (3.2)  as  a  basis  for  a  model 
of  f(xl  ),  I  =  k,  k  + 1,  k  +  2,....  We  let  fn  l  =  k,  k  + 1,  A:  +  2,...  denote  our  estimates  of 
f(x'ni),  1  =  k,k  +  l,k  +  2,....  Controls  are  defined  as  ( Nk,nk ),  (Nk+l,nk+1),  (Nk+2,nk+2),... 
and  the  dynamic  equation  for  the  state  f]  is 

fM  =  +  r*  (f,  -  r*),  /  N„r  ),  l  =  k,  k  + 1,  k  +  2,...  (3.3) 

with  initial  condition  fk  =  rf  and  where  the  equality  indicates  equality  in  distribution. 

We  define  c(N,  n)  to  be  the  computational  cost  of  carrying  out  n  iterations  of  the 
solver  applied  to  the  AP  with  a  sample  of  size  N,  where  the  terminal  state  is  characterized 
by  c(l,0)  =  0.  Tenninal  states  discussed  in  Royset  (2009)  are  defined  by  X*,  and  are 
translated  for  our  case  as 

T:={fei  \%-r*  <s} . 
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The  set  of  feasible  controls  R{%)  are  defined  as:  If 

^  eT,  then  R(d;)  =  {(1,0)};  otherwise,  R(£;)  =  ¥  x¥ .  We  seek  an  admissible  stationary 
policy  r:  j  ->¥  x¥  u{0}  with  r(f,)eR(fl)  which  minimizes  the  computational  cost 
by  evaluating  a  planning  horizon  of  feasible  controls  ( Nk,nk ),  (Nk+1,nk+l), 

(Nk+2,nk+2),... .  The  sample-size  control-problem  stated  in  Royset  (2009)  which 
accomplishes  this  task  is 


JkArf’ 


,re,ra)\=  lim  sup  E 


TjAAf,)) 


l=k 


subject  to  constraints  (3.3).  Here  E  is  the  expectation  with  respect  to  the  disturbances  on 
the  stages  k,  k  +  1,...,  s  due  to  the  truncated  normal  distribution  in  (3.3)  Finally,  we 
define  the  surrogate  sample-size  control  problem  by 


S - SSCP, (rf  ,r\re,ra)\  inf  Jkr (rf , r  ,r0,ra). 


Figure  6.  Representative  illustration  of  dynamic  program 


(3.4) 


Now,  assume  we  have  estimated  f(xkQ)  =  f  (x^  )  at  the  beginning  of  stage  k:  see 
Figure  6.  We  wish  to  choose  a  pair  ( Nk,nk )  that  will  be  computationally  efficient  so  that 

at  the  end  of  the  kth  stage  shown  in  Figure  5  we  have  progressed  toward  /*  at  a 
“reasonable”  computational  cost.  As  shown  by  the  dotted  line  in  Figure  5,  we  define  a 
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tolerance  a  to  determine  a  range  of  near-optimal  objective  values  that  are  acceptable  for 
selection  as  /* .  Further,  we  assume  that  the  computational  cost  of  the  kth  stage  is 
defined  as  c(Nk,nk )  =  Nknk .  Then,  from  (3.2),  we  see  that  selections  of  Nk  and  nk  have 
varying  effects.  A  small  Nk  increases  the  variability  in  the  truncated  distribution 
whereas  a  large  Nk  compresses  the  distribution,  reducing  the  variability.  The  number  of 
solver  iterations  also  impacts  the  cost  by  affecting  the  mean  of  the  truncated  distribution 
in  (3.2).  As  nk  increases,  the  mean  of  the  distribution  moves  closer  to/*,  but  also 
increases  the  computational  cost.  The  optimal  policy  of  S-SSCP^/y,?-*,^,^)  balances 
the  computational  cost  of  selecting  large  Nk  and  nk  and  the  likelihood  that  the  CA 

reaches  a  near-optimal  objective  value  in  the  current  stage.  Hence,  we  expect  that  policy 
to  be  reasonably  efficient  even  though  it  is  based  on  several  approximations. 

To  solve  S-SSCPi(r/,r*,r6l,rT)  approximately,  we  discretize  the  state  space  and 

the  truncated  nonnal  distribution  and  then  apply  backward  recursion  to  the  resulting 
dynamic  program.  We  refer  to  the  policy  computed  in  this  manner  as  Model-Predictive 
Control  (MPC).  We  adopt  the  discretization  technique,  as  in  Royset  (2009),  on  a 
planning  horizon  of  10  future  stages.  A  general  depiction  of  the  discretization  scheme 
applied  to  Figure  6  is  shown  in  Figures  7  and  8.  An  example  of  a  relatively  low  cost 
selection  of  Nk  and  nk  is  shown  in  Figure  7.  Here,  a  relatively  small  Nk  with  a  choice  of 

small  nk  provides  little  gains  toward  /* .  Although  progress  is  made  toward  /* , 
additional  stages  of  CA  are  necessary  arrive  at  a  near-optimal  value. 
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Figure  8  depicts  a  situation  where  a  large  Nk  and  a  large  /?/f  provide  a  near-optimal 
objective  value  close  to  f*  but  at  a  high  computational  price.  For  this  situation,  future 
stages  of  CA  may  not  be  necessary  in  order  to  achieve  a  near-optimal  objective  value. 

B.  PARAMETER  ESTIMATION 

We  rely  on  estimates  of  the  parameters  f(xkn~x),f*,  6 ,  and  cr2(x* )  in 
S-SSCPjt(r/.,r*,7"6),rT)  as  we  describe  next. 
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1.  Estimating  Variance 


Upon  completion  of  nk  iterations  of  the  solver  with  a  sample  of  size  Nk  in  stage 
k,  the  set  of  iterates  {xf}"*0  and  function  values  {fN  (xf)}"*0  are  known.  We  use  this 
information  to  estimate  parameters  rf,  r* ,  r0,  and  r  2 .  First,  we  discuss  our  estimate  for 
r  2 .  We  let  &k+l  denote  our  estimate  of  a1  (x* )  for  stage  k  +  1  and  set  it  equal  to  the 
sample  variance  at  the  last  iterate,  i.e., 


r 

<7 


2 


N„~ 


tZ(^« 
1  /= 1 


V)-A«))2 


2.  Estimation  of  Rate-of-Convergence  Coefficient  and  Optimal  Objective 
Value 

Next,  we  determine  rg,  the  estimate  of  the  rate-of-convergence  coefficient  9  of 

the  nonlinear  programming  solver  used  in  Step  3  of  CA.  We  adopt  the  method  in  Fie  and 
Polak  (1990)  to  estimate  9 ,  but  modify  it  slightly,  adding  an  exponential  smoothing  step 
to  avoid  large  changes  in  the  estimate.  As  the  estimate  of  9  is  updated  after  each  stage, 
we  let  9k  be  the  estimate  of  9  available  at  the  beginning  of  stage  k  and  set  rg  =  9k  in 
S-SSCP^rV,,^). 

As  in  Royset  (2009),  the  (biased)  estimate  of  /*  at  the  beginning  of  the  kth  stage 
is  computed  by  a  weighted  average  of  estimates  of  ,  /  =  l,2,...,k-\ ,  and  is  denoted  by 
fl .  Unlike  the  approach  used  by  Royset,  f*  is  computed  with  a  fixed,  conservative 

estimate  of  the  rate  of  convergence  denoted  by  9  .  The  meaning  of  “conservative”  is 
discussed  in  more  detail  in  the  explanation  of  Step  6  of  Subroutine  PE  below. 

Subroutine  PE (k,  9k  fk  ,  {fNk  (xf  )}*) 

Parameters:  Exponential  smoothing  parameter  t//  e  (0, 1 ) ,  conservative  estimate 
9  e  (0, 1)  of  rate  of  convergence  9  and  tolerance  sg>  0  . 
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2. 


Input  data:  Previous  stage  index  k;  estimates  9k  and  ft  ;  function  values 
from  stage  k. 

Output:  Returns  estimates  0k+1  and  ft+x . 

Step  1 :  Set  9  =  0k. 

Step  2:  Estimate  the  minimum  of  /’v(x*  )  by: 

b_  i  y ./;  (<)  ^ 

nkh  \-0"k~l 

Step  3:  Solve  least-square  problem: 

(a,b)  =  arg  min  Y  (log(/„  (x*  )  -  b)  -  i  log  a  -  bf 

«.*  i=0 

Step  4:  If  |  6 -a  \<  sg,  set  6  =  a  and  go  to  Step  5.  Else,  set  6  -  a  and  go  to  Step 


Step  5 :  Set  0k+l  =  y/0  +  (1  -  y/)9k 

Step  6:  Compute  conservative  estimate  of  the  minimum  value  of  fN(xkn) 


m,,  :=  mm 


fNS<k)-onk-lfNM) 


i=0,l,. ../it — 1 


1-0' 


Step  7:  Estimate  /* : 


:= 


Step  8:  Return  9 k  ,  and 


z;:>- 


— 7 - m.  - ‘j1 — -  ft 


From  our  assumption  of  a  uniformly  linearly  convergent  solver,  it  follows  that 
fNt  >  (fNt  (< ) - 0n^fNt  (xf)) / (1  - 0”>-‘)  for  all  i  =  0,  1,2,  ...,  nk  - 1.  Step  2  averages 
these  lower-bounding  estimates  with  the  current  estimate  of  6*  and  uses  this  as  an 
estimate  of  ft  .  With  this  estimate,  we  compute  an  estimate  of  the  rate-of-convergence 
coefficient  in  a  least-squares  sense.  That  is,  we  use  log-linear  regression  to  compute  the 
rate-of-convergence  coefficient  to  best  fit  the  sequence  {/^  (xf)}"*0.  We  define  a 

tolerance  for  computing  the  rate-of-convergence  coefficient,  denoted  by  sg ,  which 
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determines  when  to  terminate  the  regression  technique.  Next,  we  perform  exponential 
smoothing  so  as  not  to  let  the  estimate  of  9  vary  too  much  from  one  stage  to  another. 
The  value  mk  in  Step  6  of  Subroutine  PE  is  guaranteed  to  be  a  lower  bound  on  f*  only 

if  0  >9 .  Hence,  we  recommend  9  be  set  to  a  value  close  to  1 .  Since  f,*+]  is  simply  the 
weighted  average  of  m,,l  =  1,2,..., k ,  and  E[f*\<f*  (see  Chapter  II,  Section  A, 
Subsection  3),  we  therefore  find  that  f*+1  is  a  lower  bound  on  /* ,  on  average,  when 
9>9. 


C.  STOPPING  CRITERION 

We  could  implement  a  stopping  criterion  based  upon  a  hypothesis  test  of  Karush- 
Kuhn-Tucker  (KKT)  conditions  as  developed  in  Ruszczynski  and  Shapiro  (2007).  As 
stated  in  Ruszczynski  and  Shapiro  (2007),  suppose  that  the  feasible  set  X  is  defined  by 
equality  and  inequality  constraints  in  the  form 

A:={xe  j  " :  g,(x)  =  0,  i  =  g,(x)<  0,  i  =  q  +  l,...p}, 

where  the  g,  (x)  are  smooth  deterministic  functions.  If  only  equality  constraints  are 
present  and  the  gradient  vectors  Vg((x),  i  =  l,...,q  are  linearly  independent,  then  the 

hypothesis  test  of  KKT  conditions  can  be  based  upon  an  asymptotically  noncentral  chi- 
square  distribution.  If  the  assumption  of  linearly  independent  gradient  vectors  cannot  be 
met,  a  degenerate  solution  is  presented  due  to  redundancy  in  the  constraint  functions.  In 
the  case  of  both  equality  and  inequality  constraints,  a  similar  result  is  available  (see 
Ruszczynski  and  Shapiro,  2007),  which  also  relies  on  the  linear  independence  of 
gradients  of  active  constraints.  Since  such  linear-independence  assumptions  may  often 
fail  in  practical  application,  we  have  chosen  to  adopt  an  approach  using  optimality  gaps 
as  defined  by  Mak  et  al.  (1999),  and  similarly  by  Ruszczynski  and  Shapiro  (2007).  We 
base  our  stopping  criterion  upon  this  approach  with  a  small  modification. 
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With  the  iterate  jc*  found  at  the  completion  of  stage  k,  we  again  estimate  /(jc*  ) 
as  before,  with  /  (jc*  )  using  a  new  independent  sample  of  size  N* .  Here  we  elect  to 
use  a  large  sample  size  to  obtain  an  accurate  approximation  of  f(xk„k)  .  Calculation  effort 
is  not  substantially  increased  by  this  procedure  as  we  are  not  performing  an  optimization. 
From  the  central  limit  theorem,  a  probabilistic  upper  bound  on  /*  is  approximately 
normally  distributed  with  mean  /(jc*  )  and  variance  cr2(xk  )  /  N*  for  large  N* . 


The  modification  of  the  method  described  in  Mak  et  al.  (1999)  occurs  in  our 
construction  of  a  lower  bound  as  described  in  Subroutine  PE.  While  Mak  et  al.  (1999) 
use  the  average  of  the  optimal  values  of  a  set  of  APs,  we  construct  a  lower  bound  on  /* 
by  averaging  lower  bounds  on  optimal  values  of  the  APs  for  each  stage.  Our  method 
tends  to  be  more  conservative  as  it  is  based  upon  an  assumption  of  a  rate  of  convergence. 
From  Royset  (2009),  we  see  that  our  lower  bound  is  approximately  normally  distributed 

with  mean  /*  and  variance  N,,  for  large  Nv N2,..., Nk  and  nvn2,...,nk . 

Hence  the  inequality 


Prob  [f{xknt)<  f  +  s)>  O 


(3.5) 


holds  approximately.  We  therefore  stop  the  calculations  when  the  right-hand  side  in  (3.5) 
exceeds  a  selected  confidence  level  8 ,  typically  0.95  or  0.99. 
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IV.  COMPUTATIONAL  STUDY 


This  chapter  presents  a  computational  study  of  two  network-flow  problems  to  test 
how  well  Model-Predictive  Control  reduces  computation  times  compared  to  alternative 
policies  for  selecting  sample  sizes.  We  define  two  network-flow  problems  with  random 
congestion  and  present  results  of  Model-Predictive  Control  as  compared  to  heuristic 
control  of  sample  sizes. 

A.  SINGLE-COMMODITY  NETWORK  FLOW 


To  develop  a  single-commodity  flow  problem  (SCF)  for  testing,  we  consider  the 
generic  congestion  model  for  single-commodity  flows  as  described  by  Ahuja  et  al., 
(1993,  p.  651),  but  modify  it  to  include  random  congestion.  Here,  the  generic  model  has 
a  nonlinear  objective  function  of  the  form 


£ 


xu 


mm 

xgX  A/f  —  jr 

(yM  LV±  i]  AiJ 


(4.1) 


where  M..  is  the  nominal  capacity  of  arc  (/,  j)  and  xtj  is  the  flow  of  a  single  commodity 

on  arc  (/,/).  For  review  of  commodity  flow  and  congestion  modeling  we  refer  the  reader 
to  pages  109-124  in  Heam  et.  al.  (2001),  and  to  Marcotte  and  Nguyen  (1998)  and 
Bergendorf  et  al.  (1997).  We  first  present  an  SCF  problem  and  then  advance  to  a  multi- 
commodity  flow  problem.  Even  though  the  SCF  problem  is  a  special  case  of  the 
multiple-commodity  flow  problem,  we  present  SCF  first,  due  to  the  relative  ease  of 
explaining  this  simpler  problem. 


We  consider  a  graph  G  =  ( N ,  A),  where  N  and  A  are  sets  of  nodes  and  arcs, 
respectively.  The  specific  graph  considered  in  this  study  is  shown  in  Figure  9.  This  grid 
network  can  flow  commodity  left-to-right,  north-to-south  and  south-to-north,  but  not 
right-to-left.  The  test-problem  grid  has  50  nodes  and  134  arcs.  The  start  node,  denoted 
by  5,  is  the  supply  source  and  the  terminal  node,  t,  is  the  demand  sink.  Individual  costs 
associated  with  the  flow  across  an  arc  are  assigned  as  random  numbers  from  a  normal 
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distribution  with  parameters  that  will  be  specified,  and  we  define  a  congestion  parameter 
on  each  arc  by  generating  a  log-normally  distributed  random  variable,  with  parameters 
that  will  be  specified. 


Figure  9.  Transportation  Grid  Network 


We  formulate  an  SCF  problem  as  follows: 

Indices 

i,j  nodes,  i,jeN 

C iJ )  arcs  (i,  j)  e  A 

Data 

My  capacity  of  arc  (/,/)• 

C..  unit  cost  to  ship  commodity  on  arc  (/,  /). 

jUj  mean  of  log-normal  random  variable  denoting  congestion  on  arc  (/,  /). 

variance  of  log-normal  random  variable  denoting  congestion  on  arc  (/,  /). 
D  demand  of  commodity. 
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Random  Variables 


(Oy  congestion  parameter  on  arc  (/,  /);  this  is  a  log-normal  random  variable, 

with  mean  /Uy  and  variance  ov . 

Decision  Variables 

Xj  amount  of  commodity  shipped  on  arc  (/,/)• 

Mathematical  Formulation 


min 


CUXij 


y  - 

(i,j)sA  ( 1  +  O)  -  )Mj  Xy 


s-t.  I  x..-  X  A/ 

j:(j,i)eA  i:(i,j)^A 


-D  if  i  =  s 
<  0  if  i^N\{s,t} 

D  if  i  =  t 


0  <  Xy  <  My  , 


V 


hj 


(4.2) 


We  note  that  the  expectation  in  the  objective  function  can  be  computed  by 
evaluating  \A\  one-dimensional  integrals,  and  thus  a  simpler  method  for  solving  this 
model  is  available.  However,  this  model  serves  as  a  simple  example  to  illustrate  our 
solution  approach,  which  applies  to  more  general  situations. 

We  assign  500  units  of  supply  at  node  5,  with  a  corresponding  500  units  of 
demand  at  node  t.  Arc  capacities  are  chosen  as  M.  =100  for  all  (/',/) .  Based  on 

preliminary  numerical  experiments,  we  find  that  6  =  0.993  is  sufficient  to  obtain  lower 
bounds  on  f*  in  Step  6  of  Subroutine  PE. 


B.  MULTI-COMMODITY  NETWORK  FLOW 

For  testing,  we  also  consider  a  congestion  problem  for  a  multi-commodity  flow 
problem  in  a  transportation  network  (MCF),  and  reuse  the  grid  network  described  in 
Subsection  A.  The  formulation  is  as  follows: 
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Indices 

i,j  nodes,  i,jeN 

( ij )  arcs  (/,/)  e  A 

p  commodity,  p  e  {1,2,...,  P) 

Data 

M.j  capacity  of  arc  (/,/). 

C'l  unit  cost  to  ship  commodity  p  on  arc  (/,  /). 

mean  of  log-normal  random  variable  denoting  congestion  on  arc  (/,  /). 

crl  variance  of  log-normal  random  variable  denoting  congestion  on  arc  (/,  /). 

Dp  demand  of  commodity  p. 

Random  Variables 

(Oy  congestion  parameter  on  arc  (/,  /);  this  is  a  log-normal  random  variable, 

with  mean  //..  and  variance  crfj . 

Variables 

xfj  amount  of  commodity  p  shipped  on  arc  (/,/). 

Mathematical  Formulation 


min 


E 


I 


_ £ _ 

(i+ffls)Ms-y 


P 


s.t. 


-Dp  ifi  =  s 

<  0  if  i  e  N\{s,t},  V  p 
Dp  if/  =  t 


(4.3) 


ZX'7  “  MJ  ’  V  hi 
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®-xij  ’ 


v  i,j,p 


For  this  problem,  arc  capacity  is  increased  to  MtJ  =  150  for  all  (/,/)  to  allow  for 

increased  flow  from  additional  commodities.  We  consider  two  commodities  with 
supplies  at  5  equal  to  500  and  300  and  demands  at  t  equal  to  500  and  300,  respectively. 
In  MCF,  preliminary  experimentation  shows  that  9  =  0.997  tends  to  provide  a  valid 
lower  bound  of  f*  in  Step  6  of  Subroutine  PE  and  we  adopt  that  value  for  6  . 

C.  COMPUTATIONAL  STUDY 

For  this  computational  study,  we  apply  the  parameters  described  next  to  both  SCF 
and  MCF.  We  use  the  PGM  nonlinear-programming  algorithm  with  Armijo  step  size 
rule;  for  example,  see,  Polak  (1997,  p.  67)  and  Bertsekas  (1999,  p.  31).  The  quadratic 
direction- finding  problem  in  the  PGM  is  solved  using  LSSOL  (Gill  et  al.,  1986)  as 
implemented  in  TOMLAB  7.0  (Holmstrom,  2008).  We  use  parameters 
a  =  0.5  and  ft  =  0.8  in  the  Armijo  step-size  rule  and  in  Subroutine  PE  use  an  exponential 
smoothing  parameter  y/  =  1  /  3  and  tolerance  sg  =  0.0001 . 

For  stopping  criterion,  we  draw  a  new  independent  sample  of  size  N *  =  10000  to 
evaluate  fN,  (xknf )  and  use  a  stopping  confidence  level  of  S  =  0.95 .  We  use 

(jUy,<jy)  =  ( 3,4)  for  all  (/,/)  as  parameters  for  the  log-nonnal  distributed  random  variable 
(Oy  representing  congestion.  Arc  costs  C..  and  C’’  are  generated  from  a  normal 
distribution  of  random  numbers  with  mean  80  and  standard  deviation  20.  Additionally, 
we  set  the  relative  optimality  tolerance  to  0.01  for  use  in  calculations,  i.e.,  s  =  0.01  /’.( , 
on  stage  k. 

For  comparison  studies,  we  consider  two  versions  of  Model-Predictive  Control; 
MPC1  and  MPC2.  In  MPC1,  Model-Predictive  Control  is  applied  to  all  stages  of  the 
conceptual  algorithm,  including  the  first  stage.  We  find  empirically  that  MPC  might 
have  poor  control  in  the  initial  stage  when  the  parameters  estimated  in 
S-SSCPi(r/.,r*,r6l,rT)  are  inferior  estimates.  Hence,  we  also  consider  MPC2,  where 
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Model-Predictive  Control  is  used  from  the  second  stage  of  the  conceptual  algorithm.  The 
first  stage  uses  a  predetermined  sample  size  and  number  of  solver  iterations.  We 
examine  three  choices  for  the  first-stage  policy  resulting  in  the  following  three  cases  of 
MPC2: 

1 .  MPC2a.  nx  =450,  Nx  =  100 ,  and  remaining  stages  use  MPC. 

2.  MPC2b.  nx  =  600,  /V,  =  100 ,  and  remaining  stages  use  MPC. 

3.  MPC2c.  n}  =  900,  /V,  =100,  and  remaining  stages  use  MPC. 

Choices  of  nx  for  these  cases  of  MPC2  are  detennined  by  solving  6'h  =  0.01 ,  0"'  =  0.05, 

and  6 =0.1,  where  6  is  the  conservative  rate-of-convergence  coefficient  as  discussed 
in  Chapter  3,  Section  B,  Subsection  2. 

As  a  basis  for  comparison,  we  consider  the  following  heuristic  policies: 

1 .  Fixed  policy.  Predetermine  Nk  and  nk  and  keep  fixed  throughout  each 
stage  of  CA. 

2.  Additive  policy.  Predetermine  nk  and  add  a  predetermined  number  to  Nk 
at  the  beginning  of  each  stage  of  CA  (i.e.,  Nk+l  =50  +  Nk). 

3.  Multiplicative  policy.  As  in  additive  policy,  predetennine  nk,  and  adjust 
sample  size  by  a  multiplicative  factor  at  the  beginning  of  each  stage 
(i.e.,  Nk+X  =l.2Nk) . 

We  use  an  initial  sample  size  Nx  =10  for  all  heuristic  policies,  except  for  the 

fixed  policy  for  which  Nk  =0.5  •  N*  =5000  for  all  k.  In  order  to  use  the  PGM  within 

CA,  we  must  first  find  an  initial  feasible  solution  for  both  SCF  and  MCF  to  start  the 
calculations.  We  do  so  by  formulating  and  solving  the  following  linear  program  in  the 
case  of  SCF: 
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SCF-LP:  min  V  C  x 

J  i] 

(i,j)£A 

-D  if i  =  s 

s.t.  ^  Xjf-  ^  xtj  =  <  0  if  i  e  A  \  {s,  t} 

1  W  [  /)  jf/-/ 

0  <  Xy  <  My  ,  V  i,7 

and  the  following  linear  program  in  the  case  of  MCF: 

MCF-LP:  min  Z  ZCW 

-Dp  ifi  =  s 

s.t.  z  xJi-  Z  xv=<  0  Vp 

-I  [pf  if/=/ 

Z-Fy  -  Mu  ’  v  hi 

P 

0<xp  ,  V  i,j,p 

We  implement  our  network-flow  problems  in  Matlab  Version  7.7.0  on  a  desktop 
computer  running  Windows  XP  with  3.73  GHz  processor  speed  and  3.25  GB  of  RAM. 
SCF-LP  and  MCF-LP  are  solved  to  find  an  initial  feasible  solution  for  both  SCF  and 
MCF  using  the  linear  programming  solver  linprog  in  the  optimization  toolbox. 

For  comparison  studies,  we  record  the  computing  time  of  CA  using  each  of  the 
MPC  policies  with  the  computing  time  for  each  of  the  other  policies  considered.  We 
evaluate  each  of  the  heuristic  policies  with  a  different  predetermined  control  on  the 
number  of  solver  iterations.  Evaluations  are  run  with  nk  for  all  k  set  at  5,  25,  50,  75  and 

100  iterations.  In  the  additive  approach,  the  sample  size  is  increased  by  100  at  the 
beginning  of  every  stage.  For  the  multiplicative  approach,  two  separate  heuristics  are 
considered.  First,  we  evaluate  the  policy  with  an  adjustment  to  sample  size  as 
Nk+x=\.5Nk  for  all  k  and  then  increase  the  adjustment  control  on  sample  size  to 

Nk+l  =  2 Nk  for  all  k. 
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The  results  summarized  in  Table  1  provide  average  computing  times  over  20  runs 
of  the  CA  with  standard  deviations  for  the  SCF  problem.  The  first  column  lists  the 
individual  policies  mentioned  above  for  determining  (Nk ,  nk ) .  The  second  and  third 

columns  give  the  average  and  standard  deviations,  respectively,  of  the  total 
computational  times  to  reach  a  near-optimal  objective  value  in  the  SCF  problem. 


Policy 

SCF  Computational  Times  (sec.) 

avg  st  dev 

MPC1 

23.17 

3.60 

MPC2a,  «!  =  450, N x=  100 

14.87 

2.67 

MPC2b,  n  !  =  600,  =  100 

18.38 

2.48 

MPC2c,  n  !  =  900,  N  x  =  100 

24.77 

3.15 

Fixed,  n  =  5 

443.30 

2.48 

Fixed,  n  =25 

638.12 

42.84 

Fixed,  n  =50 

605.40 

50.24 

Fixed,  n  =75 

614.86 

38.27 

Fixed,  n  =100 

631.84 

50.60 

Additive,  n  =  5 

398.07 

16.78 

Additive,  n  =  25 

87.58 

7.78 

Additive,  n  =  50 

49.26 

6.76 

Additive,  n  =75 

41.54 

6.81 

Additive,  n  =100 

40.81 

9.08 

Mult  1.5,  n  =5 

>  1100 

Mult  1.5,  n  =25 

81.95 

13.54 

Mult  1.5,  n  =  50 

31.20 

6.31 

Mult  1.5,  n  =  75 

29.85 

7.85 

Mult  1.5,  n  =  100 

27.95 

5.93 

Mult  2.0,  n  =  5 

>  1100 

Mult  2.0,  n  =  25 

>  1100 

Mult  2.0,  n  =  50 

77.06 

7.52 

Mult  2.0,  n  =  75 

35.98 

7.34 

Mult  2.0,  n  =100 

32.56 

8.29 

Table  1.  Average  and  standard  deviation  of  computing  times  of  CA  with  Model-Predictive 
Control  (MPC)  policies  and  heuristic  policies  applied  to  SCF. 


In  Table  1,  the  best  heuristic  policy  with  respect  to  computational  time  is  the 


multiplicative  policy  with  a  multiplicative  factor  of  1.5  and  ^  =100  for  all  k.  In 

comparison,  MPC1  finds  a  near-optimal  objective  value  nearly  five  seconds  faster  and 
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does  so  with  a  39.3%  reduction  in  standard  deviation  of  computational  time  over  the  20 
independent  runs.  MPC2a  improves  further  still,  offering  a  reduction  in  computing  time 
of  nearly  47%.  Additionally,  the  standard  deviation  between  the  independent  runs  drops 
55%  as  compared  to  the  best  heuristic  policy.  The  computational  times  recorded  do  not 
reflect  the  time  required  to  determine  the  Model-Predictive  Control,  i.e.,  to  solve 
approximately  S-SSCPt(r/,r*,/'0,r£7) .  We  elect  to  exclude  this  time  because  for  large, 

real-world  problems,  computing  times  for  the  minimization  calculations  and  checking 
stopping  criterion  are  expected  to  be  considerably  larger  than  computing  times  for 
determining  (Nk ,  nk ) . 

Several  policies  considered  for  SCF  return  results  that  are  costly  regarding 
computing  times.  In  those  cases,  we  terminate  the  calculations  after  1100  seconds  and  do 
not  compute  averages:  see  rows  12,  17,  and  18  of  Table  1.  For  the  policies  with  times 
greater  than  1100  seconds,  the  relatively  small  number  of  solver  iterations  is  not 
sufficient  to  make  substantial  gains  towards  /* .  In  these  cases,  Nk  grows  quite  large 
and  computing  time  suffers  from  the  large  sample  size.  For  each  of  the  problems,  Nk  is 

limited  to  400,000  to  avoid  exhausting  computer  memory,  and  in  each  of  these  cases,  the 
sample  size  grows  to  this  limit. 

The  results  summarized  in  Table  2  for  MCF,  provide  average  computing  times 
over  20  runs  of  the  CA,  along  with  standard  deviations. 
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Policy 

MCF  Computational  Times  (sec.) 

avg  std  dev 

MPC1 

40.70 

7.21 

MPC2a,  «!  =  450,  N  x=  100 

26.46 

4.22 

MPC2b,  n  i  =  600,  N  i  =  100 

28.84 

4.35 

MPC2c,  n  !  =  900,  N  x  =  100 

37.57 

2.95 

Fixed,  n  =  5 

>  1000 

Fixed,  n  =  25 

891.65 

47.96 

Fixed,  n  =50 

905.67 

61.53 

Fixed,  n  =75 

905.11 

58.92 

Fixed,  n  =100 

891.89 

57.60 

Additive,  n  =  5 

612.40 

41.05 

Additive,  n  =25 

143.12 

14.81 

Additive,  n  =  50 

79.86 

9.91 

Additive,  n  =75 

66.74 

10.75 

Additive,  n  =100 

54.86 

7.24 

Mult  1.5,  n  =  5 

>  1000 

Mult  1.5,  n  =  25 

151.96 

33.66 

Mult  1.5,  n  =  50 

48.66 

8.86 

Mult  1.5,  n  =  75 

44.41 

7.51 

Mult  1.5,  n  =  100 

48.12 

12.71 

Mult  2.0,  n  =  5 

>  1000 

Mult  2.0,  n  =  25 

>  1000 

Mult  2.0,  n  =  50 

142.55 

46.12 

Mult  2.0,  n  =  75 

61.95 

16.52 

Mult  2.0,  n  =  100 

52.64 

11.67 

Table  2.  Average  and  standard  deviation  of  computing  times  of  CA  with  Model-Predictive 
Control  (MPC)  policies  and  heuristic  policies  applied  to  MCF. 


As  in  SCF,  a  number  of  the  policies  make  the  sample  size  grow  until  it  hits  the 
limit  of  400,000,  thereby  affecting  overall  computing  times.  In  these  cases  we  terminate 
the  calculations  after  1000  seconds  and  do  not  compute  averages. 

The  best  heuristic  policy  for  MCF  is  again  a  multiplicative  policy.  However  in 
this  larger  problem,  computational  time  is  best  when  nk  =75  for  all  k :  compare  this  to 

SCF,  where  computational  time  is  best  when  nk  =  100  for  all  k.  MPC1  improves  on  this 
computational  time  by  nearly  four  seconds  and  does  so  with  essentially  the  same 

variability  of  computational  time  between  runs  as  the  best  heuristic  policy. 
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We  see  that  modifying  MPC  in  the  first  stage  gives  further  computational  savings. 
Policy  MPC2a  shows  an  improvement  of  40%  in  overall  computing  time,  on  average,  and 
improves  the  standard  deviation  of  computing  time  by  almost  44%  over  the  20 
independent  runs.  These  results  indicate  that  while  the  MPC  typically  provides  a  “good” 
policy  for  selecting  (Nk,  nk ) ,  the  estimates  of  parameters  in  S  -  SSC  Pt  ( rf ,  r" ,  re ,  ra )  for 
the  first  stage  are  rather  poor  and  a  heuristic  policy  may  be  better  in  that  stage. 

To  verify  that  the  stopping  test  (3.5)  does  not  cause  premature  termination  of  CA, 
we  compute  a  lower  bound  on  /*  as  described  in  Mak  et.  ah,  (1999).  Specifically,  we 
run  the  PGM  on  the  AP  with  N  =  10000  until  that  algorithm  stalls  and  record  the  last 
function  value.  This  is  an  estimate  of  f*N .  We  repeat  this  process  30  times.  By  the 
central  limit  theorem,  the  average  of  these  function  values  is  approximately  normal  and 
provides  a  lower  bound  /  on  /* .  We  find  that  in  all  160  runs  of  the  MPC  policies,  the 
probability  that  the  last  solution  found  is  no  worse  than  (1  +  0.01)  /  is  essentially  1.0. 

Hence,  the  stopping  test  (3.5)  is  rather  conservative,  as  zero  unsatisfactory  solutions  is 
well  within  the  0.05  •  160  =  8  expected  when  8  =  0.95 . 

MPC1  solves  each  of  the  network-flow  problems  faster  than  any  of  the  heuristic 
policies  considered.  Additional  reductions  in  computing  time  are  gained  with  MPC2 
where  the  first  stage  of  the  CA  is  primarily  used  to  estimate  parameters  in 
S-SSCP^.rV,,^)  fork >2. 
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V.  CONCLUSIONS 


This  thesis  develops  an  efficient  optimization  algorithm  for  approximately  solving 
stochastic  nonlinear  programming  problems  whose  objective  functions  are  sample- 
average  approximations.  We  demonstrate  improvement  in  computation  times  by 
approximately  solving  a  discrete-time  optimal-control  problem  to  select  a  policy  of  well- 
balanced  sample  sizes  and  number  of  solver  iterations  for  each  stage  of  the  algorithm. 
This  policy,  referred  to  as  the  Model-Predictive  Control  policy  (MPC),  is  compared 
against  alternative  heuristic  policies  for  selecting  sample  sizes  and  solver  iterations. 
MPC  approximately  solves  a  single-commodity  network-flow  problem  up  to  17%  faster, 
on  average,  than  the  best  heuristic  policy.  Furthermore,  the  optimal-control  problem 
provides  a  40%  reduction  in  standard  deviation  of  computing  times  over  a  set  of 
independent  runs  of  the  algorithm  on  identical  problem  instances.  When  we  fix  the 
number  of  solver  iterations  in  the  first  stage  and  then  proceed  with  MPC,  we  improve  the 
computing  time,  on  average,  by  nearly  47%  and  reduce  the  standard  deviation  between 
runs  by  more  than  one  half. 

The  application  of  the  discrete-time  optimal-control  problem  to  a  larger  multi- 
commodity  network-flow  problem  shows  an  8.4%  improvement,  on  average,  in 
computational  time  over  the  best  heuristic  policy  with  essentially  the  same  variation  of 
overall  computational  time  between  the  20  independent  runs.  With  the  first-stage 
modification  to  Model-Predictive  Control,  we  improve  computing  time  by  40%,  on 
average,  compared  to  the  same  heuristic  policy  and  reduce  standard  deviation  between 
runs  by  44%. 

The  algorithm  developed  to  solve  nonlinear  stochastic  programs  shows 
considerable  promise  and  offers  significant  potential  for  further  study. 
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